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Abstract 

This study employed entropy-based second law post-processing analysis to characterize the various 
thermodynamic losses inside a 3-space solution domain (gas spring + heat exchanger + regenerator) operating 
under conditions of oscillating pressure and oscillating flow. The 3-space solution domain is adapted from the 
2-space solution domain (gas spring + heat exchanger) in Kornhauser’s MIT test rig by modifying the heat 
exchanger space to include a porous regenerator system. A thermal non-equilibrium model which assumes 
that the regenerator porous matrix and gas average temperatures can differ by several degrees at a given 
axial location and time during the cycle is employed. An important and primary objective of this study is the 
development and application of a thermodynamic loss post-processor to characterize the major 
thermodynamic losses inside the 3-space model. It is anticipated that the experience gained from 
thermodynamic loss analysis of the simple 3-space model can be extrapolated to more complex systems like 
the Stirling engine. It is hoped that successful development of loss post-processors will facilitate the 
improvement of the optimization capability of Stirling engine analysis codes through better understanding of 
the heat transfer and power losses. It is also anticipated that the incorporation of a successful thermal non- 
equilibrium model of the regenerator in Stirling engine CFD analysis codes, will improve our ability to 
accurately model Stirling regenerators relative to current multi-dimensional thermal-equilibrium porous- 
media models. 
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I. 


Introduction 


A. Background 

Under the auspices of NASA's Radioisotope Power System (RPS) development program, multiple efforts are 
currently underway, both in-house at NASA Glenn Research Center (GRC) and under various grants and contracts, 
to develop a high-efficiency Stirling Radioisotope Generator (SRG) for possible use on future NASA Space Science 
Missions. The SRG is being developed for multi-mission use including providing electric power for unmanned Mars 
rovers and deep space missions in environments with and without atmospheres. One of the multiple efforts 
underway includes the development of a multi-dimensional Stirling computational fluid dynamics code, including 
second law analysis post-processing to separate various thermodynamic losses 34,37 in Stirling engines. The second 
law loss analysis effort documented in this paper supports the efforts underway in the Stirling research community 
and is sponsored by NASA grant NNC05AA24A. 

Heat engines convert heat to useful work. This conversion process is inherently irreversible due to the non-ideal 
nature of power systems. Internal system irreversibilities resulting from viscous friction, heat transfer and non- 
equilibrium processes destroy part of the available work of a power system resulting in thermodynamic losses which 
impact engine thermodynamic performance adversely. There are two major design objectives for space-power 
applications of reciprocating heat engines like the Stirling engine: (1) To maximize performance by minimizing 
thermodynamic losses — in order to minimize fuel requirements (expensive radioisotopes for example), and (2) To 
minimize system mass — in order to minimize propulsion fuel requirements. In order to improve the engine 
thermodynamic performance, it is necessary to identify and minimize the features guilty of the available energy loss 
within the system. Thermo-fluid system loss analysis and management is clearly an area of study that is generating a 
great deal of interest. 

Optimal engine performance requires good heat transfer to and from the working fluid and maximum 
conversion of the input heat to useful work by reducing thermodynamic losses as noted above. Until recently the 
heat transfer estimates used in reciprocating heat engines have been based on extrapolation of expressions that were 
developed for steady-pressure, steady-flow conditions. Because the ordinary, steady-state convective model contains 
no term to account for the oscillatory effect in variable volumes, designs of reciprocating engines are likely not 
based on these steady flow expressions. Using two Massachusetts Institute of Technology (MIT) test rigs (gas spring 
and gas spring+heat exchanger) Kornhauser 20 confirmed experimentally that there was a fundamental difference 
between steady and oscillatory flow heat transfer in the variable cylinder volumes of the MIT test rigs. Building on 
the work done by Lee, Smith, Faulkner, and Chafe 6 ’ 13 ’ 21 , Kornhauser started the development of expressions suitable 
for oscillating pressure and oscillating flow conditions. He closely integrated experiment with analysis to achieve 
useful results with good qualitative but limited quantitative success. 

In this study, entropy-based second law post-processing analysis is employed to characterize the various 
thermodynamic losses inside a 3 -space solution domain (gas spring + heat exchanger + regenerator) operating under 
conditions of oscillating pressure and oscillating flow. The 3 -space solution domain was adapted from the 2-space 
solution domain (gas spring+heat exchanger) in Kornhauser’ s test rig 20 by modifying the heat exchanger space to 
include a porous regenerator system. A thermal non-equilibrium model which assumes that the regenerator porous 
matrix and gas average temperatures can differ by several degrees at a given axial location and time during the cycle 
is employed. A survey of the porous-media literature supports the need for thermal non- equilibrium porous-media 
models for Stirling regenerators ’ ’ ’ . Numerical simulation results of the pressure, temperature, velocity and 
surface heat transfer variations in the 3 -space domain are post processed for thermodynamic loss calculations. These 
simulation and loss calculation results are compared with similar results of an earlier study 11,12 in order to evaluate 
the effect of adding a regenerator to the 2-space model. Also results of surface heat flux, temperature difference 
between the gas temperature at the radial center of the heat exchanger and the heat exchanger wall, temperature 
contours and velocity vectors are compared to relevant results from the literature. 

Second-law analysis focusing on entropy generation has been playing a dominant role in recent times and has 
proven to be a useful tool in identifying the mechanisms and system components that are responsible for 
thermodynamic losses and for indicating how to minimize these losses in practical equipment in order to improve 
performance 10 . Entropy generation destroys part of the available work of a system and is associated with 
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thermodynamic irreversibilities related to pressure drop, finite heat transfer, friction, mixing and other nonidealities 
within systems. Past studies have described various analytical and empirical techniques for entropy-based 
optimization of engineering systems. Bejan 2,4,5 has focused on the different mechanisms behind entropy generation 
in applied thermal engineering. Numerous other investigations (mostly numerical) have been performed to 
determine entropy generation and irreversibility profiles for different geometric configurations, flow situations, and 
thermal boundary conditions 24,25,28 . 

A valuable tool in the design of a high performance engine is a numerical code with multi-dimensional 
modeling capability and an ability to closely simulate the thermal-fluid processes inside the engine and account for 
all the thermodynamic losses. This study utilized 1-D Sage and 2-D Fluent commercial numerical codes to obtain 
computer models of the 3-space solution domain and to analyze and post-process the numerical simulation results. A 
thermodynamic loss post-processor is developed to characterize the major thermodynamic losses inside the 3 -space 
model components. 

Sage is a 1-D, multi-variable thermodynamic modeling package that supports simulation and optimization of 
spring-mass-damper resonant systems and Stirling-cycle coolers and engines. In the Sage code, an engineering 
system is built up from component parts obtained from component palettes. The components function as a whole by 
virtue of their connections which could be due to gas flow, heat flow, pressure, density, etc. Sage calculations are 
performed via its solver and optimizer routines. 

Fluent is a state-of-the-art computational fluid dynamics (CFD) software package with comprehensive and 
flexible physical modeling and multi-physics capabilities for simulating fluid flow and heat transfer in complex 
geometries. The Fluent package includes the solver (FLUENT), the preprocessor (GAMBIT) for geometry modeling 
and mesh generation, an additional preprocessor (TGrid) that can generate volume meshes from existing boundary 
meshes and filters (translators) for import of surface and volume meshes from a variety of CAD/CAE packages. All 
functions (specifying problem type and numerical technique, setting boundary and initial conditions, defining fluid 
properties, etc.) are accessible in FLUENT through the Graphical User Interface (GUI) an interactive, menu-driven 
interface with many tools to compute a solution, visualize the flow physics, animate transient data sets, as well as to 
extract data for post processing the numerical results. Tools such as point, line, rake and surface are used to obtain 
the data of required parameters. 

B. Scope/Objectives 

The rest of this paper can be broadly categorized as: 

• An Overview of the Stirling Engine 

• Modeling of the 3 -space Solution Domain 

• Regenerator-Effect of change from Sage 2-space to Sage 3-space Modeling 

• Thermodynamic Loss Modeling 

• Numerical Simulation Results 

• Thermodynamic Loss Post-Processing 

• Conclusions and Suggestions for Future Work 

In order to fully understand the fundamental objective of this research effort, it is necessary to provide an 
overview of the Stirling engine and its operation. A description of Komhauser’s MIT 2-space Test Rig is next 
presented so as to provide the proper context for the 3 -space solution domain used this study. Following this, the 1-D 
Sage and 2-D Fluent computer modeling procedures for the 3-space solution domain are presented. After that, the 
effect of adding a regenerator to the 2-space solution domain in Kornhauser’s MIT test rig is discussed. Next, 
conservation equations of mass, momentum and energy are united with the second law of thermodynamics to derive 
a complete set of entropy generation equations appropriate for thermodynamic loss modeling. The numerical 
simulation results of the fluid pressure, temperature, fluid velocity and surface heat transfer variations are presented 
and post-processed for thermodynamic losses. Finally, conclusions are drawn and suggestions for future work 
provided. 

An important and primary objective of this study is the use of the second law analysis for the development and 
application of a thermodynamic loss post-processor to characterize the major thermodynamic losses inside the 3- 
space model. It is anticipated that the experience gained from the thermodynamic loss analysis of the simple 3 -space 
model can be extrapolated to more complex systems like the Stirling engine via the incorporation of loss post- 
processors in Stirling engine analysis codes in order to improve its optimization capability. It is also anticipated that 
a thermal non-equilibrium model of the regenerator such as that presented here, when incorporated in Stirling engine 
analysis codes, will improve our ability to accurately model Stirling regenerators relative to current thermal- 
equilibrium porous-media models. 


NAS A/CR— 2008-2 1 5480 


3 



II. 


An Overview of the Stirling Engine 


Figure 1 below shows a schematic representation of the Stirling Technology Demonstration Convertor (TDC), a 
55-watt space power Stirling engine prototype developed by Infinia Corporation (formerly Stirling Technology 
Corporation (STC)). It is essentially a free piston machine that generates electric power from a thermal energy input. 
The convertor is of the gamma type, gas-coupled, and single-acting. It can be divided into two basic subsystems - 
electro-mechanical and thermo-fluid. 



Figure 1. Schematic of STC’s Stirling Convertor [Courtesy: www.nasa.gov l. 


The electro-mechanical subsystem consists of a pressure vessel containing a flexurally supported power 
piston/linear alternator assembly. The power piston separates the thermodynamic working space in the thermo-fluid 
subsystem to the right of the piston from the “bounce” space, a weak but measurable gas spring between the piston 
and the casing. 

The thermo-fluid subsystem can be considered to have two distinct fluid circuits (internal and external) for most 
applications. The internal engine thermodynamic circuit is filled with working fluid (usually a gas: air, helium or 
hydrogen) at an elevated pressure and is comprised of a flexurally supported displacer/rod assembly, two variable 
volumes (compression and expansion volumes) and three heat exchangers in series: the heater, the regenerator and 
the cooler. The heater-head channels the heat from an external heat source, an electric heater (radioisotopes are the 
planned heat source for space application), into the expansion space and heats the working fluid to a temperature of 
about 650°C. The regenerator is a duct packed with some porous matrix (solid containing interconnected pores for 
fluid flow). It is often made by stacking together a number of fine wire mesh screens, random fibers (for the TDC), 
foam metal plugs, or perforated disks or felts to form a kind of metallic heat sponge. The cooler, in contact with a 
heat sink, a water circuit in the laboratory (solid conduction to a radiator for space power applications), extracts 
excess heat from the compression space providing a cold-sink temperature of about 80°C. 

The external circuit (not shown in Figure 1) provides thermal energy which could come from just about any 
thermal energy source, available at a sufficiently high temperature. Stirling engines have been run on concentrated 
solar energy, thermal energy storage batteries, metal combustion, isotope nuclear energy, as well as a variety of 
liquid and solid fuels (oil, coal, gas, etc.). 

Clearance seals used for the power piston and displacer isolate the gas working spaces. The radially- stiff 
flexures make tight clearance seals possible. The convertor uses hot-end materials capable of operating at the hot- 
end temperature over the planned life of the machine. 

The TDC operates in a closed regenerative thermodynamic cycle created in the thermo-fluid subsystem. During 
each working cycle, the power piston and the displacer reciprocate in a coordinated almost sinusoidal fashion (with 
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the power piston motion lagging the displacer motion 16 ), shuttling the gas at different temperature levels through the 
heat exchangers between the expansion and compression spaces. Pressures acting on each end of the displacer are 
only slightly different, relative to the mean pressure level. On the other hand, the relatively large “bounce” space to 
the left of the power piston (see Figure 1), compared to the working space volume (compression, expansion and heat 
exchanger volumes, etc), ensures that the pressure excursions in the bounce space will be substantially less than 
inside the engine working space. The larger pressure excursions inside the working spaces help drive the piston, 
alternately compressing and expanding the gas. The power piston’s attached alternator parts oscillate within the 
alternator's magnetic field to produce electric power. The power piston extracts net power from the working space 
gas due to the larger pressure excursions on that side. The oscillating gas pressures within the working spaces, the 
spring-mass-damper systems of the piston and displacer, and the oscillating load on the piston produced by 
interaction with the alternator, all interact to produce sustained oscillating motion of the piston and displacer — when 
engine geometries, piston and displacer masses and heat sources and sinks are appropriately chosen. 

The action of the Stirling engine regenerator is most easily understood by first imagining the gas behavior 
without a regenerator. In the absence of the regenerator, hot gas would be transferred by the action of the displacer 
directly from the expansion space/heater into the cooler/compression space, where it would have to be cooled. The 
heat extracted during the cooling would be rejected and lost. When the gas is subsequently returned to the expansion 
space, it would have to be reheated drawing more heat from the heat source. Extra heat would therefore be added 
and rejected during the cycle, with a consequence of loss of efficiency. With the regenerator in place, the hot gas 
gets cooled gradually as it passes from the hot end to the cold end, by giving up the heat to the regenerator solid 
matrix and therefore the gas leaves the regenerator already cooled, minimizing the heat to be rejected in the cold 
space. On the return journey, the gas is gradually heated up as it moves up the temperature gradient towards the 
expansion space, by picking up heat from the solid matrix that was deposited during the previous cycle. Thus the gas 
emerges into the heater with considerable heat already added, minimizing the heat to be added by the external 
source 24 . The regenerator thus serves as an economizer for storing heat during one part of the engine cycle and for 
re-use during another part. The inclusion of an effective regenerator results in a substantial increase in Stirling 
engine efficiency and is necessary for the Stirling engine to be of practical use. 

Extensive efforts have been focused on improvement in regenerator technology. These efforts have been 
categorized into areas of materials and geometry, numerical modeling, and experimental measurement ’ ’ ’ ’ . A 
NASA regenerator research grant effort led by Cleveland State University, with subcontractor assistance from the 
University of Minnesota, Gedeon Associates, and Sunpower Inc. has been providing computational and 
experimental results to support definition of various empirical parameters and “closure” relations needed in defining 
a thermal, non-equilibrium, macroscopic, porous-media model for use in multi-D Stirling codes for regenerator 
simulation 33 . 
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III. Modeling of the 3-space Solution Domain 

(Komhauser’s MIT 2-space Solution Domain Modified). 


A. Kornhauser’s MIT 2-space Test Rig 



Table 1 . Gas Spring + Heat Exchanger Dimensions 20 


Cylinder Bore (Diameter) 

50.80 mm (2.0 in) 

Piston Stroke 

76.2 mm (3.0 in) 

Volume Ratio 

2.0 

Annulus Outside Diameter 

44.5 mm (1.75 in) 

Annulus Inside Diameter 

39.4 mm (1.55 in) 

Annulus Wall-to-Wall Distance 

2.5 mm (0.10 in) 

Annulus Eength 

445 mm (17.5 in) 

Min. Cyl.-Head to Piston Clearance 2.9 mm (0.11 in) 


Figure 2. Komhauser’s MIT Cylinder + HXer Test Rig. 

Figure 2 shows a schematic representation of Kornhauser’s piston-cylinder+heat exchanger test rig. A piston- 
cylinder device is mounted on a compressor base. The compressor piston drives the piston in the 50.8 mm (2.00 in.) 
diameter cylinder. The piston is sealed to the cylinder with a buna-n O-ring located more than a stroke’s length from 
the piston face so that frictional heating of the cylinder wall would only affect results minimally. The piston top 
surface is brass, while the cylinder wall and head are made of steel. The fixed cylinder head has an annular opening 
leading into an annular dead-ended heat exchanger space such that flow can continuously pass between the heat 
exchanger and cylinder as the cylinder piston expands and compresses the gas similar to that in a Stirling engine 
compression space volume. The test section consists of the gas spring + heat exchanger space. A piston stroke length 
of 76.2 mm (3.00 in.) and volume ratio of 2.0 are used in this study. A volume ratio is defined to be the maximum 
cylinder volume (piston at “bottom-dead-center” (BDC)) over the minimum cylinder volume (piston at “top-dead- 
center” (TDC)). 

The heat exchanger annulus is 44.5 mm (1.75 in) outside diameter and 39.4 mm (1.55 in) inside diameter, for a 
wall-to-wall distance of 2.5 mm (0.10 in). The annulus is 445 mm (17.5 in) long, so that a volume ratio of 2.0 for the 
combined cylinder and heat exchanger resulted in a very small cylinder clearance volume; the cylinder-head 
distance is nominally 2.9 mm (0.11 in) at top center position. The heat exchanger entrance has the same cross- 
sectional dimensions as the heat exchanger itself, the entrance comers being as sharp as could easily be machined. 
The inner wall of the heat exchanger space is made of steel. The outer wall of the heat exchanger space is steel lined 
with a 1.98 mm (0.078 in) layer of bronze-filled TFE (Teflon). This lining was chosen for thermal properties 
matching those of the Pyrex glass substrate of the surface temperature transducers. Because of the extended heat 
transfer surface in the annular heat exchanger, the energy flows are more complicated than in a simple gas spring. 
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The apparatus is belt driven by a D.C. motor to allow for speed adjustment. The apparatus fill line is a capillary 
tube of negligible volume. The gas spring + heat exchanger dimensions are tabulated in Table 1. 

The 3-space model solution domain (gas spring + heat exchanger + regenerator) was adapted from the 2-space 
model solution domain (gas spring+heat exchanger) for Komhauser’s 2-space test rig (Figure 2) by modifying the 
heat exchanger space to include a porous regenerator matrix. 1-D Sage and 2-D Fluent commercial numerical codes 
are used to model, analyze and post-process the thermal-fluid phenomena inside the 3 -space domain. The modeling 
procedures using Sage and Fluent are presented in subsections III.A.l and III.A.2 below. Note that in Figures 3 
(Sage model) and 4 (Fluent model), the regenerator is shown located at the heat exchanger end opposite from the 
end adjacent to the cylinder. The justification for this is discussed in Section IV. Note also that this 3-space device is 
different from the usual Stirling engine. The regenerator is not located between a heater and a cooler as in a Stirling 
engine. Also, whereas the 3 -space device contains a single variable volume on one side of a heat-exchanger- 
regenerator arrangement, the Stirling engine contains a variable volume on each side of a cooler-regenerator-heater 
heat-exchanger circuit. The flow dynamics in the 3-space solution domain is expected to approximate some of the 
flow dynamics in the Stirling engine thermo-fluid subsystem 

1. Sage 3 -space Model (1-D) 

Figure 3 is a representation of the Sage graphical interface illustrating the 1-D model of the 3-space solution 
domain. 



Figure 3. Sage Model of the 3-space Solution Domain. 


Sage provides 8 menus - File, Display, Edit, Scan, Specify, Process, Options, and Help - each with self- 
explanatory commands and graphical edit windows with model-component palettes. When “show window” in the 
edit menu is activated, the highest level graphical edit window (Stirling machine window) and its model-component 
palettes appear. Each model-component palette contains buttons arranged in tab-selected pages corresponding to 
different categories of model components. Figure 3 shows the Stirling machine window with six Stirling engine 
model component palettes - “Basic”, “Canisters”, “Heat Exchanger", “Phasr Moving Parts”, “Gt Moving Parts” and 
“Composite”. A Sage model is, in part, a collection of component parts assembled in a graphical edit window by 
clicking and dragging the components into the window and connected in a particular way to form a complete 
system. Model components are organized in a hierarchical structure. The root model component contains a number 
of sub-components which may themselves contain sub-sub-components. Child (sub) models are accessed by double- 
clicking on parent models. The labeled arrows sticking out of the sides of components are boundary connectors. 
They are joined together as matched pairs with numbers (2, 3, 4, etc.) identifying the match. The labels are meant to 
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suggest what information is transferred across the connection between components. Sage is also able to perform 
such functions as open and save files, increase/decrease connector level, copy, cut and paste components. 

The 3 -space solution domain (gas spring + heat exchanger + regenerator) is modeled using Sage’s Stirling 
machine model component palettes. The variable volume cylinder space in Fig. 2 and the charge pressure inside it 
are modeled in the Sage code using a generic cylinder (renamed “gas spring cylinder” in Fig. 3) and a pressure 
source, both obtained from the “Basic” component palette at the Stirling machine level. The pressure source comes 
with a built in steady-state density connection (p s td y ) and acts as an isobaric gas reservoir in that the density in the 
gas domain adapts itself so that the time average pressure is continuous across the connection. The cylinder-space 
gas and isothermal surface, child models of the gas spring cylinder, are obtained from the “Gas Domain” and 
“Cylinder Walls” component palettes respectively at the gas spring cylinder level. The arrows labeled 2 in the gas 
spring edit window indicate a space-time varying heat flow connection Q Gxt between the z-face of the cylinder-space 
gas and the isothermal surface. The connection arrows from the positive-facing volume displacement phasor, gas 
charge line and negative gas inlet, all child models of the cylinder-space gas, are moved up to the Stirling machine 
level for pressure connections (P p h S r) to the constrained piston (label 3), density connections (p st d y ) to the pressure 
source (label 4) and mass flow connection to the heat exchanger (label 7) respectively. The positive-facing volume 
displacement phasor represents the volume change of the gas space due to the motion of the piston. The gas charge 
line and negative gas inlet are obtained from the “Charge/Inlef ’ component palette in the cylinder- space gas level. 

The piston in Figure 2 is modeled using the constrained piston obtained from the “Phsr Moving Parts” 
component palette at the Stirling machine level. Its child model, the negative facing-area phasor, is obtained from 
the “Mechanical Attachment” component palette in the constrained piston level. Its pressure connection arrow is 
moved up to the Stirling machine level for connection there to the gas spring cylinder (label 3). 

The annular canister in Figure 3 is used to model the heat exchanger and the regenerator. It is obtained from the 
“Canisters” component palette in the Stirling machine level. The generic matrix, child model of the heat exchanger, 
represents the internal structure of the heat exchanger. It is obtained from the “Matrices” component palette in the 
annular canister level. The matrix gas and isothermal surface, child models of the generic matrix, are obtained from 
the “Gas Domain” and “Matrix Solids” component palettes respectively in the generic matrix level. The arrows 
labeled 2 in the generic matrix edit form indicate a heat flow connection between the matrix gas and the isothermal 
surface. The negative and positive gas inlets, child models of the matrix gas, are obtained from the “Charge/Inlef ’ 
component palette in the matrix gas level. The mass flow connection arrows from the negative and positive gas 
inlets are moved up to the Stirling machine level for connection there to the regenerator (label 8) and gas spring 
cylinder (label 7). The regenerator is modeled in a similar manner as for the heat exchanger except for the 
requirement of a rigorous surface child model for the surface condition of the random fiber matrix representing the 
internal structure of the regenerator. The positive gas inlet, child model of the matrix gas, models the gas flowing 
from the regenerator into the heat exchanger (label 8). We see from Figure 3 and the foregoing descriptions that the 
3-space model components communicate with each other using density (p st d y ), mass flow (m Gt ), heat flow (Q Gx t) 
and pressure (P p h S r) boundary connections. 

Once a model structure is created in Sage, numerical inputs for the model component are either specified or 
modified. Sometimes user-defined variables (special output variables) can be added to model components. Sage 
numerical simulation is initiated after initializing component and overall model parameters. Sage calculations are 
performed via its solver and optimizer routines. 

2. Fluent 3-space Model (2-D) 

The Fluent 2-D model of the 3 -space solution domain is illustrated in Figure 4 below. The origin (0,0) of the x-y 
coordinate system is located at the intersection of the cylinder head and the line of symmetry. The Fluent code 
exploits the symmetry of the problem domain to model only one-half of the domain. The geometry of the physical 
domain of interest is represented within the preprocessor GAMBIT using fundamental geometric entities (points, 
lines, arcs, circles, curves, splines and surfaces) which can be manipulated (translated, rotated, scaled, projected, 
split, joined, etc.) as desired. Grid generation, the process of discretizing the problem domain with individual cells 
over which the flow equations are integrated, follows the geometric representation of the domain of interest. The 
locations of the comer points of these cells constitute the “grid” or “mesh” which is stored in a data base. Fluent 
provides complete mesh flexibility, including stmctured and unstmctured meshes, different mesh types and ability to 
refine and coarsen the grid based on the flow solution. The power law grid distribution type is used to obtain a finer 
grid at the boundaries of the model for more accurate resolution of the flow features. The baseline grid size 100x20 
is presented in Figure 4 for illustration. The optimum linear dimension of the regenerator was determined via a 
systematic parameter analysis using the Sage code. This is discussed in Section IV. 
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After modeling is completed in GAMBIT, the meshed geometry is exported to Fluent where the problem is 
completely formulated using an interactive, menu-driven Graphical User Interface (GUI). Problem formulation 
involves specifying the problem type, model options, volume conditions, boundary conditions, initial conditions, 
formulation type (implicit or explicit), the numerical technique needed to solve the flow equations and solution 
controls. The instantaneous position X p (t) and velocity U p (t) of the piston face in the x-direction are defined as 


Xp(t) = 



1 + sin 



( 1 ) 


U P (t) 


dX p (t) 

dt 


S 

— cocos 
2 


r 


V 



( 2 ) 


where parameter L c is the length of the cylinder clearance volume when the piston is at TDC, co is angular frequency, 
t is time and S is the piston stroke. Fluent uses the velocity equation. A user defined function (UDF) is written to 
input the velocity equation into Fluent. After problem formulation, simulation is performed using the FLUENT flow 
solver module. 
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IV. Regenerator-Effect on Change from Sage “2-space” to “3-space” Modeling 


The results of a prior study 12 of a Sage model of the MIT 2-space solution domain (Figure 5(a)) indicated a 
cylinder-cooling effect, with heat being drawn from the environment through the cylinder walls into the cylinder 
volume and pumped from the heat exchanger gas volume through the heat exchanger cylinder walls into the 
environment. Since the 3 -space solution domain used in this study was adapted from the MIT 2-space solution 
domain by modifying the heat exchanger space to include a porous regenerator system, it was decided to investigate 
the regenerator impact on the cylinder-cooling effect in changing from Sage 2-space to 3 -space modeling. Table I 
summarizes the results of the investigation. 

The second row of Table 1 shows results obtained using a Sage model of the original MIT 2-space solution 
domain. This model assumed isothermal wall temperatures of 294 K for both the heat exchanger and the cylinder. 
Since the purpose of a “real cooler” is to reduce the temperature of the cooled area relative to the environment, the 
temperature of the cylinder walls was reduced (with everything else remaining the same) over a range of 
temperatures below 294 K, in order to obtain the wall temperature below which no cooling results. The effects of 
changing the cylinder wall temperature are shown in Runs 1 thorough 5 (Table 1). Based on these results, the 
cylinder temperature (289 K) for the next series of runs was established in Run 3. Run 3 (also designated Run A) 
was thus chosen to be the reference “cooler” case for subsequent runs B through H. 

Next, in preparation for converting part of the 44.5 cm long heat exchanger into a regenerator, the heat 
exchanger was first split arbitrarily into two parts (40 cm part and 4.5 cm part) such that of the 16 equilength cells in 
the 44.5 cm long heat exchanger, 14 equilength cells are now contained in the 40 cm part and 2 equilength cells in 
the 4.5 cm part sectioned out from the end of the heat exchanger facing the cylinder for use as the regenerator. With 
no regenerator porous matrix added, Run B was then initiated with ensuing results essentially the same as for Run A. 
The slight differences are likely due to the non identical equilength cells. 

The 4.5 cm part of the heat exchanger was next converted into a regenerator resulting essentially in a 3-space 
Sage model of the solution domain - cylinder, regenerator and heat exchanger. The regenerator matrix is modeled as 
either a generic random-fiber or a “dedicated” random-fiber. The regenerator wall is adiabatic and heat exchanger 
walls is kept fixed at 294 K. The regenerator position relative to the heat exchanger and cylinder is now varied for 
best possible location based on cooling performance. With the regenerator between the heat exchanger and the 
cylinder (Figure 5(b); Runs C and D, Table 1) the cooling effect of the device was completely “wiped out”, 
apparently because the regenerator acted as a “heat dam” to prevent heat from moving from the cylinder to the heat 
exchanger. Runs C and D differ only in the representations of the 90% porosity regenerator matrix (generic-matrix- 
based random- fiber model in Run C, and “dedicated” random- fiber model for Run D). 

In Runs E and F (90% porosity generic-matrix-based random fiber model and dedicated 90% porosity random- 
fiber model, respectively), the regenerator was placed at the end of the heat exchanger away from the cylinder 
(Figure 5(c)). A significant improvement in cooling performance, relative to the performance of the original 2-space 
Sage model (Run A) and Run B was noted. The coefficient of performance (COP) improved from -0.145 to - 0.21 1. 
It is also observed that the switch from generic matrix based random fiber to “dedicated” random fiber matrix (Run 
C to Run D and Run E to Run F) has less effect with the regenerator on the heat exchanger end away from the 
cylinder (Runs E and F) than with the regenerator between the heat exchanger and the cylinder (Runs C and D) - — 
perhaps because “less is going on” at the heat exchanger end away from the cylinder. 

In Run G, the impact of the 90% porosity matrix at the end of the heat exchanger opposite the cylinder was 
checked by removing the generic matrix based 90% porosity random fiber matrix. Very similar results to Runs A and 
B were obtained, as anticipated. 

In Run H, it was thought necessary to investigate if the improvements shown in Runs E and F - regenerator at 
end of heat exchanger away from cylinder — were due primarily to the 10% reduction in the volume produced by 
adding the 90% porosity random fiber. This investigation was done by simply cutting off 4.5 cm of the heat 
exchanger. Removing this 4.5 cm’s “worth of volume” does improve the cooling performance relative to Runs A and 
B (from COP -0.146 to COP - 0.153), but not nearly as much as putting a 90% porosity random fiber matrix in the 
4.5 cm. part of the heat exchanger away from the cylinder (COP - 0.146 to COP - 0.211). Thus the impact of the 
regenerator matrix on the cooling performance is significant but it is not presently clear why the addition of the 
regenerator matrix at this end improves the performance. Since the regenerator is not located between a heater and a 
cooler, as in a Stirling engine, it is not as obvious how storage and release of heat to and from the regenerator matrix 
is improving the cooling efficiency. 
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Table 1: Summary Results of Regenerator-Effect (i.e. Results of Changing from Sage 2-space to 3-space Modeling) 


Comments 

Sage model of MIT’s 2-space 
test rig. (Figure 5(a)) 

Next 5 mns: Reduce Cylinder Temperature 



Run A (Reference) 



Rest of mns involve splitting-off part of heat exchanger - - to serve as regenerator, or not 

Serves as a reference case for runs B 
through H 

Essentially same result as Run A. Slight 
difference is likely due to different 
equilength cell distribs. 

Regenerator acts like a “heat dam”. 
Most of heat goes out through the 
cylinder wall (Fig. 5(b)) 

Almost same results as for Run C. Slight 
differences due to differences in random 
fiber matrix model (Fig. 5(b)). 

The switch from generic to dedicated 
random fiber matrix has less effect with 
the regenerator at end of HXer away 
from the cylinder than with regenerator 
between HXer and cylinder (Fig. 5(c)). 

Dramatic improvement in cooling 
performance - significant when 
compared to Runs A and B (Fig. 5(c)). 

Results very similar to Runs A and B. 

Removing this 4.5 cm’s “worth of 
volume” does improve cooling 
performance relative to Runs A ana 
B, but not nearly as much as for 
90% porosity random fiber matrix 
in the 4.5 cm. 
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Figure 5(b): Sage Model Schematic Corresponding to Runs C and D 
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Figure 5(c): Sage Model Schematic Corresponding to Runs E and F 
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Figure 6: COP vs. Regenerator Length (m) 
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With the 3 -space model configured as in Runs E and F, a series of parameter optimization studies can be 
designed and executed in an effort to obtain optimum values for regenerator parameters such as regenerator length, 
wire diameter, annulus diameter (inner and outer) and porosity with the objective being the maximization of the 
cylinder-cooling effect. The Run F configuration of the 3 -space model is used in this study because of its slight COP 
advantage over the Run E configuration (see Table 1). The results of coefficient of performance (COP) values over a 
range of regenerator lengths, which appears to be one of the easiest optimizations to do, are presented in graphical 
form in Figure 6. The variation in the regenerator length is made by increasing its length from 4.5 cm. while 
ensuring that the overall length of the regenerator-heat exchanger domain remains constant at 44.5 cm. The optimum 
cylinder-cooling effect (COP = 0.33019) is observed at the regenerator length of 22 cm. Other regenerator parameter 
optimization studies will be left for some future considerations. 

In order to observe the effect of the regenerator on the numerical simulation results of the original 2-space 
model 12 , the 3-space model operating conditions were set to match those of the 2-space model (201.7 RPM, 1.008 
MPa T wa ii - 294 K). For this set-up, the best possible cylinder-cooling effect (COP - 0.69647) is observed at a 
regenerator length of 20 cm. as illustrated in Figure 6. The optimum length of the regenerator is thus set at 20 cm. 
for the remaining numerical simulations to be conducted with the cylinder wall temperature of 294 K. A discussion 
of the effect of the regenerator on the numerical simulation results is presented in Section V. 
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IV. Thermodynamic Loss Modeling 

The fluid flow in the 3 -space domain is highly compressible due mainly to internal volume changes. The 
competing effects of the heat transfer due to temperature gradients and flow friction points to thermodynamic loss 
models which require the coupling effects of the mass, momentum and energy conservation equations be united with 
the concept of entropy generation which relies on the second law of thermodynamics. 

A. Local or Microscopic Equations 

The theoretical development of the thermodynamic loss models begins with a set of generalized local or 
microscopic conserved equations of mass, momentum and energy and non-conserved entropy equation for the two 
mechanically distinct phases, a (flowing fluid) and (3 (stationary solid) inside the non-porous regions of the 3 -space 
domain. It is assumed that the a and (3 phases do not react chemically; that the fluid phase is assumed to behave as a 
linearly viscous fluid (Newtonian) whose properties are known functions of temperature and pressure, or other 
combinations of state properties; that the density of the solid is constant and that the no-slip assumption at the fluid- 
solid interface is valid. The following local or microscopic equations of mass, momentum, energy and entropy that 
govern the flow and describe the losses in the non-porous regions of the 3 -space domain are similar to those derived 
in any fluid mechanics text. 

1. Conservation of Mass Equation 

The conservation of mass equation, applicable only to the a-phase, is: 


dpq 

a 


+ V '(pq U q) = 0 


(3) 


This equation is a scalar equation with 4 unknowns: the gas density p a and three components of the gas velocity U a . 


2. Conservation of Momentum Equations 


The conservation of momentum equation is also applicable only to the a-phase. The differential statement of the 
conservation of momentum for a Newtonian fluid with no body forces can be written as: 


d(pqUq) 

dt 


V '(p« U a U a) + V P 



f-v(nqV-u a ) -(v-p a v)u a 
2(vn a -v)u a - (v n a x Q a ) = 0 


(4) 


The parameter Q a is the vorticity vector. This equation is a vector equation with three components corresponding 

to the three components of velocity. The conservation of momentum equation produces 1 additional unknown, the 
fluid pressure, p a . The fluid dynamic viscosity, p a , is a material property. For many fluids, p a depends significantly 
on temperature and when appreciable temperature differences exist in the flow field it is necessary to regard p a as a 
function of position. 

3. Conservation of Energy Equation 

The energy equation is applicable to both the a and (3 phases. 


d(pqUq) 

0t 


a-phase 

V-(p a u a h a -k a VT a )= 0 


(5a) 


p-phase 

<5b) 


The energy equations introduce 3 additional unknowns; the fluid and solid temperatures (T a and Tp), and the 
fluid specific internal energy u a - The thermal conductivities k a and kp, like the viscous coefficient p a , are 

temperature dependent material properties. The parameters cp, and (pc)p, are the solid specific heat and heat capacity 
per unit volume respectively. 

The gas enthalpy, h a = u a + (p a /p a )• Although not immediately obvious from Eq. (5a), it can be shown 26 that 
the gradient of the enthalpy flow term, V-(p a u a h a ) = V (p a u a u a ) + p a V u a - O a where the viscous 
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dissipation function O a is defined as: 


O a 


— (2/ 3) p. (V • u) 2 + 2[jJ 



+ 


8v^ 

dy. 





' 8v + 5u^ 2 
v 5x + 5y y 


( dw 

dvY 

' dn 

dw^i 

— 

+ — + 

— 

+ — 


dz ) \ 

K dz 

3x J 


( 6 ) 


4. Equations of State 

Thus far, the system of six equations (mass, three components of momentum and the two phases of energy) 
containing eight unknowns (p a , three components of u a , p a , T a , T p , u a ) requires two additional equations to close 
the problem. These equations are the thermodynamic equations of state which can be expressed generally as 


Pa = /l(P.»“«) (7) 

T„=/2(p„,ft„) (8) 

The particular equations of state are arbitrary, in the sense that the particular forms do not affect the form of the 
governing equations. Ideal or real gas equations can be used. 

5. Entropy Generation Equation 

The differential statement of the second law of thermodynamics 3 can be re-written in different forms for the a 
and (3 phases as: 


a-phase 


p-phase 


s 


fff 

gen, a 


d(pa S a) 

dt 


+ V- 

% 

+ V '(pa S a U a) ^ 0 (9a) 




T 


gen,P ^ 

l T P y 


s 


m 

gen, a 


a -VT O 

a _|_ a s, q 

T 2 T 

A a a 


(9b) 


* gen, p 


q P -vTp 

T 2 

A P 


> 0 


(9c) 


(9d) 


The entropy generation rate for the non-porous regions of the 3 -space domain is thus: 

s'" = s'" + s'" >0 nm 

gen, sys (non-porous) a, gen p,gen v lv v 

The second law postulates the existence of entropy, S, a non-conserved thermodynamic property of state that 
can be created via a generation or production term, s"' en . The parameter q represents the heat flux vector. The 

inequality indicates that the entropy generation is always positive except for totally reversible processes, in which 
case, it is zero. Equations (9) and (10) clearly show that entropy generation is due to system irreversibilities resulting 
from heat transfer, mass flow and viscous friction. 

The foregoing equations represent the generalized, non-volume averaged, microscale equations for a single 
fluid phase flowing and interacting with a stationary solid phase. 
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B. Volume-Averaged Equations 

The flow geometry in the porous regenerator is too complicated to allow for a direct application of the local or 
microscopic equations. The regenerator presents the need for volume-averaging 38 , an analytical tool for describing 
the phenomena of compressible flow and heat transfer in a porous medium. The local or microscopic conservation 
equations for the solid and fluid phases are volume averaged over a representative elementary volume (see Fig. 7 
below) to arrive at the corresponding “macroscopic” transport equations. 



Figure 7. A two-phase model of a porous medium. 

The a-phase is a Newtonian fluid and the (3-base is the solid part of the porous matrix 18 . 


The parameter n a is the unit normal vector on the a-(3 interface pointing into the (3-phase and x is the position 
vector. The subscript on the position vector indicates the phase at that position. The porous matrix is assumed to be 
stationary and rigid with known thermal properties which are functions of temperature. 

The details of the volume averaging method are discussed in Ref. 18. The resulting simplified equations are 
shown below: 


1. Volume- Averaged Conservation of Mass Equation 


f(p„r + v.(<p„)*<g„}«)=o 


(ii) 


Note that this simplified macroscopic conserved equation of mass is similar in form to the corresponding 
microscopic equation (Eq.(3)). 


2. Volume- Averaged Conservation of Momentum Equations 



I 

'V a + n a -vV-u a 

+ ^ a VU a 


V 3 ) 



+ e a 


:(%«)“ -v)u a )“ + (%«)“ 


x/Q„ 




= 0 


( 12 ) 


The parameter s a = V a /V is the gas phase volume fraction, or porosity. V a is the gas phase volume and V is the 
arbitrary averaging volume containing the porous medium (see Figure 7 above). The notations ( ) and A indicate a 


volume average and a spatial deviation quantity respectively. 
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3 . Volume-Averaged Conservation of Energy Equation 
a-phase 


(3-phase 


^ ({p. )"{“.>”)+ V • (<P, )" (u„>* (h, )■ - 

N,k„V{T„) - ) - »,h((T,)' - (T.)‘) = 0 

(13a) 


(p c )i 


5V 

dt 


l k P x P V ( T P 


+ — a„H' 


(( T p)'-( T a)“)=0 


P 


(13b) 


Volume averaging in this case results in the film heat transfer term, a term additional to the standard form of the a- 
and (3-phase microscopic energy equations (5a, 5b). Gedeon 14 recognized N k as axial conductivity enhancement due 
to thermal dispersion having a Peclet dependence of the form (ci + c 2 Re m Pr n ) where Ci and c 2 are constants and m 
and n are positive constants. This functional form for the dispersion agrees with other empirical predictions 32,34 . The 
parameters Ov and H are the surface area per unit volume and the convective heat transfer coefficient respectively. 


4. Volume- Aver aged Entropy Generation Equation 


a-phase 



(14a) 

After considerable simplification using the Maxwell 
equations the entropy generation equation reduces to 18 : 


(3-phase 



(14b) 

relations for a pure substance 3 , momentum and energy 


a-phase 



(15a) 


(3-phase 

Ae", 8 e„ >“ = k p 

V 

8 a a v H 

e e< T «> a < T A 


v ( t p) p 

J 

|< T P> P - ( T a>“ 


(15b) 


Thus the entropy generation rate equation for the gas-matrix or porous system can be written as 

(Cn,sys (porous)) = £ a ( S a'ge„ ) + £ P ( S p"ge„ ) * 0 (16) 

Substituting Eqs. (15a and 15b) into Eq. (16) the entropy generation equation for the porous regenerator system can 
be written as 


( S g'en,sys)= £ a N k k a 


V ( T J 


\ 2 


+ epkp 


f i \bA 2 

v/t r 


8„a„H 


V X y V \ p / J 

Gas conduction Matrix conduction 


( T a)“(\ 


(T p ) P -(T a )“f-s a ^fv(p a )“ > 0 


(17) 


Film heat trans fer viscous and inertial 
losses 


Eq. (17) indicates that the sources of irreversibility are gas conduction, matrix conduction, film heat transfer and 
viscous and inertial losses. 
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The volume-averaging technique assumes that the velocity, pressure and temperature can be represented in 
terms of a single large-scale averaged quantity in regions having significantly different mechanical and thermal 
properties. This assumption provides the basis for the porous-media models available in two CFD codes, Fluent and 
CFD-ACE, used at GRC for modeling Infinia’s Stirling TDC. The standard models assume that the porous matrix 
and the fluid are in local thermal equilibrium at each spatial interface between them. This is believed to be a poor 
assumption for the oscillating-flow environment within Stirling regenerators 29,34 . 1-D regenerator models, used in 
Stirling engine design, use thermal non-equilibrium models and suggest regenerator porous matrix and gas average 
temperatures can differ by several degrees at a given axial location and time during the cycle 33 . A survey of the 
porous-media literature supports the need for thermal non-equilibrium porous-media models for thermal storage 
applications such as Stirling regenerators ’ ’ ’ ’ . It is anticipated that a thermal non- equilibrium model such as that 
presented here, when incorporated in commercial CFD codes, will improve our ability to accurately model Stirling 
regenerators relative to current thermal-equilibrium porous-media models. 


C. Second Law Analysis via Entropy Generation 

Equations (10) and (17) contain several source terms which represent thermal dispersion, gas-to-matrix heat 
transfer and viscous and inertial friction. These fundamental mechanisms which are associated with internal system 
irreversibilities generate entropy within thermo-fluid systems. The second law analysis via entropy generation is 
concerned mainly with the computation of the available energy losses due to these system irreversibilities. 

Entropy generation may be calculated in one of two ways: by entropy transfer accompanied by heat transfer and 
mass flow across the external boundaries of a system (external entropy generation) or by entropy generated by 
internal processes due to heat transfer and viscous and inertial losses (internal entropy generation). In principle, the 
two methods of accounting should give the same answer. Discrepancies which usually arise are often attributable to 
numerical errors (finite-difference truncation errors, round-off errors, interpolation errors, etc.). 

1 . External Entropy Generation 

The external entropy calculation is made by integrating the heat transfer and mass flow over the surface 
between the internal calculation domain and the environment. In a reciprocating device like the MIT test rig, the 
external entropy generation over each period is obtained via a cyclic integral. 

1.1 N on-Porous Regions of the 3-space Domain: 

Substituting equations (9a) and (9c) in Eq. (10) gives 
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Expressing Eq. (18) in integral form we get 
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Taking the cyclic integral of Eq. (18) over one period gives: 
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The integral terms account for the net entropy transfer with heat and mass out of the non-porous regions of the 3- 
space domain into the surrounding environment during each cycle of operation. 


1,2 Porous Regions of the 3-space Domain: 

Substituting equations (14a and 14b) into Eq. (16) gives: 



Expressing Eq. (23) in integral form we get: 



Taking the cyclic integral of Eq. (25) over one period gives: 




The integral terms account for the net entropy transfer with heat and mass out of the porous regions of the 3 -space 
domain into the surrounding environment during each cycle of operation. 
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2. Internal Entropy Generation 

Internal entropy generation can be accounted for by tallying up the individual entropy generations due to, heat 
transfer, viscous and inertial losses in all internal processes. 

2.1 N on-Porous Regions of the 3-space Domain: 

Substituting equations (9b) and (9d) in Eq. (10) gives 


gen, sys 


qq-vT a + <k _ q P - VT P 


> o 


Expressing Eq.(29) in integral form we get 
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(31) 


Taking the cyclic integral of Eq. (31) over one period gives: 
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(34) 


Eq.(32) accounts for the internal entropy generation in a reciprocating device like the MIT test rig. From Eqs. (33 
and 34), we see that the first term on the right-hand side of Eq. (32) which involves the squares of the temperature 
gradients (see Eq. (33)) is the entropy generation due to conductive heat flow while the second term which involves 
the squares of the velocity gradients (see Eq. (34)) is the contribution due to viscous dissipation. Equation (32) 
seems to imply that the internal thermodynamic losses can be minimized by minimizing the magnitude of these local 
temperature and velocity gradients and the sum of their squares, throughout the domain of interest. Although a 
reduction in the fluid viscosity would lead to a reduction of the velocity gradients at the wall and the boundary layer 
viscous losses and an increase the fluid thermal conductivity would reduce the temperature gradients at the wall, and 
the conductive heat transfer loss; these measures are difficult or impractical to implement. It is however practical to 
minimize local velocity and temperature gradients via changes in the geometry and operating conditions such as 
frequency of operation. 
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2.2 Porous Regions of the 3-space Domain: 

Expressing the system entropy generation equation, Eq.(17), in integral form gives 
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Taking the cyclic integral of Eq.(35) over one period gives: 
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3. Availability Energy Loss: 

It has been noted above that internal entropy generation due to system irreversibilities destroys part of the 
available work of a system. The destruction of available energy impacts engine thermodynamic performance 
adversely. Thus, in order to improve the thermodynamic performance of an engineering system it is necessary to 
study the various mechanisms responsible for entropy generation in order to minimize it in engineering equipment. 

All internal entropy generations can be characterized in terms of availability energy loss defined as 3,7 : 

Availabiliy energy loss = T 0 S gen (37) 

The reference temperature T 0 is the lowest naturally occurring temperature in the system. Evaluation of the 
external entropy generation equations (22 and 28) and internal entropy generation equations (32 and 36) enable the 
corresponding availability energy loss calculations using Equation (37). Increase in entropy gives a measure of the 
extent to which the energy of a system is lost or unavailable for work during a certain process. The direct 
proportionality between the available energy loss and entropy generation provides a useful gauge to the significance 
of the various irreversibilities within engineering systems. 

The availability-loss concept allows us to think about entropy generation in terms of the more concrete notion 
of lost mechanical work. A loss in availability equates to a decrease of PV power in an engine. For example, in 
turbo-machines that generate shaft power (turbines) or absorb power (pumps, compressors), the rate of power lost 
owing to irreversibilities is proportional to a loss in availability and thus to an increase in the rate of entropy 
generation. 

4. Second Law Analysis Summary Results: 

The entropy generation results can be summarized as follows: 


External Entropy Generation: 
Non-porous domain: 
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Porous domain: 
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Internal Entropy Generation: 
Non-porous domain: 
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Porous domain: 
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Availability Energy Loss: Availabiliy energy loss = T L S gen (44) 

Loss analysis using entropy-generation rates due to heat transfer, viscous friction and non- equilibrium processes 
is a relatively new technique for assessing component performance. It offers a deep insight into the flow phenomena, 
allows a more exact calculation of losses than is possible with traditional means involving the application of loss 
correlations and provides an effective tool for improving performance. Entropy generation maps of cumulative 
amounts of all losses computed locally in the flow domain can be produced and designers can use them to detect 
critical areas (locations in which entropy generation is higher than its integral average value over the entire flow 
field). The design emphasis would then be aimed at reducing the local values for entropy production in these critical 
areas by modifying design variables so as to improve the overall performance. 

Our understanding of loss mechanisms is however far from complete. Although numerical predictions are 
valuable in predicting the heat transfer and flow structure, there are difficulties in predicting the loss accurately. This 
is due to errors in predicting the boundary layers, and laminar/turbulent transition, as well as due to false entropy 
generation due to numerical dissipation. This work provides a point of reference for incorporation of loss post- 
processors into Stirling engine numerical codes. The incorporation of a loss post-processor in Stirling engine 
numerical codes, it is believed, will facilitate the optimization of Stirling engine performance. 
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V. 


Numerical Simulation Results 


In a prior study 12 , extensive numerical simulations of the fluid flow and heat transfer phenomena under 
conditions of oscillating pressure and oscillating fluid flow inside the original MIT 2-space solution domain were 
performed using 1-D Sage and 2-D CFD-ACE+ numerical codes. In this study, 1-D Sage and 2-D Fluent numerical 
codes are used to obtain numerical simulation results in the Run F configuration (see Table 1) of the 3-space solution 
domain. As discussed in Section III, the optimum length of the regenerator is shown in Figure 6 to be 20 cm. (with a 
corresponding 24.5 cm heat exchanger length). For 2-D numerical simulations in both prior and current studies, the 
grid size, numbers of cycles of piston motion and time step were optimized prior to performing the simulations, as 
detailed in the prior study 12 . The thermal non-equilibrium porous-media models discussed in Section IV can be 
employed for the numerical simulations of the fluid flow and heat transfer effects in the regenerator. However, some 
difficulties were encountered incorporating the derived models in Section IV in the Fluent code. Following a series 
of consultations with Fluent Inc., a viable alternative for performing thermal non- equilibrium analysis was 
suggested. This entails using the TGrid feature in Fluent to generate a volume mesh of the solid phase and separating 
it from existing boundary mesh of the regenerator domain resulting in two volume meshes - a solid volume mesh 
and a fluid volume mesh. The film heat transfer term in Eq. (17) is then employed as the source term to link the two 
phases. This suggested approach was adopted in this study. 

Code (ID vs. 2D) and model (2-space vs. 3 -space) simulation results of pressure-volume diagrams and crank 
angle dependent temperature, pressure and wall heat transfer profiles are compared. To ensure a common basis for 
comparison, the reference point for the piston motion is taken at BDC (0 = 0°) moving toward TDC. Also, results of 
heat exchanger surface heat fluxes, temperature differences between the gas temperature at the radial center of the 
heat exchanger and the heat exchanger wall, domain temperature contours and velocity vectors generated by the 2D 
codes are compared with similar results obtained from the literature ’ ’ . For each comparison, important 
differences and similarities and the effect(s) of the regenerator are noted and discussed. 

1. Code (ID vs. 2D) and Model (2-space vs. 3-space) Comparisons 

/. 1 Pressure- Volume Diagrams 


2-space model 3 -space model 



Figure 8. Pressure- Volume Diagrams over a Cycle in the Cylinder Space 
(Operating Conds.: 201.7 RPM, 1.008 MPa., T w = 294 K, 480 tspc, Grid size: 147 x 46). 

Pressure-volume diagrams over a cycle in the variable volume cylinder space of the 2-space and 3-space 
domains are illustrated in Figs. 8(a,b). The indicated operating pressure is the arithmetic mean (of the compression 
and expansion space pressures) cycle pressure. Wall and piston-face temperature was fixed at 294 K. The pressure- 
volume diagrams appear qualitatively similar. Note that the maximum cylinder volume (160.4 cm 3 ), the minimum 
cylinder volume (5.9 cm 3 ) and the swept volume (154.5 cm 3 ) are the same in both models and that the maximum 
pressure in each cylinder space is on the average about 2.3 times the minimum pressure. Calculations of maximum 
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pressures at minimum compression volume compare better (~ 0.6% error) in the 3-space model than in the 2-space 
model (~ 4.3% error). On the other hand, the minimum pressures at maximum expansion volume compare slightly 
better for the 2-space model (~ 0.3% error) than for the 3-space model (~ 0.4% error). On the average, Sage and 
Fluent pressure results appear to compare better than Sage and CFD-ACE+ pressure results. 

The areas enclosed by the P-V diagrams are each equivalent to a net piston work input because the work done 
by the system on the piston during the expansion process is less than the work done on the system by the piston 
during the compression part of the cycle. The discrepancy between ID and 2D calculations of the piston net work 
input in each model was significant (~ 44% error). The ID Sage code under predicts the piston net work input. The 
inclusion of the regenerator resulted in a reduction in the minimum and maximum pressure and work-input values in 
the cylinder with the exception of the increase in Sage maximum pressure. 

1.2 Crank Angle dependent Temperature and Pressure Profiles 

2-D instantaneous temperature and pressure data are obtained at specified points (stationary and/or moving) in 
the domain of interest over a cycle of piston motion defined using crank angle, 0 , as: 

Crank Angle, 9 = [(cot + 3^r/2)- 180/^r + 270]/360 

The parameter co is angular frequency and t is time. The use of crank angle enables comparison at the same 
piston position. 

Sage employs the equation 

3 

0. = Q + ^ A n cos( ncot i + cp), i- 1,2,..., #tspc (46) 

n = 1 

to calculate the instantaneous values, ® , of the data of interest, ® , (temperature or pressure) over the number of 
time steps per cycle (#tspc). Each of the instantaneous values ® . is equal to the sum of the mean value, ® , and 
three Fourier cosine harmonics of the piston motion at each time step. The parameter A is the amplitude of ® and 
cp is the phase angle. ® ; A and cp , which are spatially averaged values over the domain are obtained from Sage 
output results. 

In the cylinder space, temperature and pressure data are obtained at specified stationary points (x,y) and moving 
points (x-node, y-node) and in the heat exchanger and regenerator spaces, temperature and pressure data are 
obtained at specified stationary points only. Note that x, y-node numbers are used as coordinates to tag the moving 
point location. The y-coordinate (0.0132673 m) or y-node (9), for example, is fixed in each domain. 

The cylinder space temperature profiles for three stationary points and three moving grid points are illustrated in 
Figs. 9(a-d) below. The x-coordinates (or x-nodes) of the points are as shown in the Figures. 

The cylinder space temperature profiles are observed to be dependent on point location and on whether the 
point is stationary or moving. Profile irregularity is more pronounced for stationary points than for moving points 
and more especially for stationary points close to the midpoint of the cylinder clearance volume (x = 0.001374 m) 
and the piston top center position (x = 0.002747 m). The sudden dip in temperature at 0 = 180° for the stationary 
point at x = 0.002747 m. may possibly be due to flow disturbances close to the piston top center position. All CFD- 
ACE+/Fluent generated peak temperatures for stationary and moving points are below Sage’s peak temperature. 
Sage over predicts the highest temperatures recorded by the 2-D codes in the 2-space and 3-space models by ~ 16K 
and ~ 11 K respectively. In general, temperature values are observed to increase with crank angle during the 
compression phase and decrease during the expansion phase of the cycle as expected. Close to the cylinder head or 
entrance to the heat exchanger (x = 0.000153 m, x-node 2), the fluid in the cylinder space experiences an 
appreciable drop in the peak temperature since almost all the fluid is pushed into the heat exchanger at this point. 
Addition of the regenerator elevates the cylinder space temperatures at stationary points especially at points close to 
the midpoint of the cylinder clearance volume (x = 0.001374 m) and piston top center position (x = 0.002747 m). 
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2-space model 


3-space model 



Figures 9(a-d). Temperature Profiles at Stationary and Moving Points in Cylinder Space 
(Operating Conds.: 201.7 RPM, 1.008 MPa, T wa n = 294 K, #tspc = 480, Grid size = 147 x 46 ). 

Temperature profiles for the stationary points in the heat exchanger space are illustrated in Figs. 9(e,f). Five and 
three stationary points are specified for the 2-space and 3 -space models respectively. Only three stationary points are 
specified in the heat exchanger space of the 3-space model because of its reduced length due to the incorporation of 
the regenerator. Profile irregularity is pronounced during the compression part of the cycle for points at or close to 
the heat exchanger entrance (0.0 and -0.11 m). As expected, maximum temperatures are recorded near the end of the 
compression process and minimum temperatures recorded close to the end of the expansion process. Unlike in the 
cylinder space, peak temperatures in the heat exchanger recorded by the 2-D codes are greater than Sage’s peak 
temperature (~ 304K @ 0 ~ 148°) and in this case, all minimum temperatures recorded by the 2-D codes are lower 
than or about equal to Sage’s minimum temperature (~ 286 K @ 0 - 274°). For the 2-space model, the highest 
temperature (~ 344. IK @ 0- 146.3°) and lowest temperature (-260. 8K @ 0 - 287.3°) are recorded near the end of 
the heat exchanger (x = - 0.33 m) disregarding the anomaly near the heat exchanger entrance. For the 3-space model, 
the highest temperature (- 333 K @ 0- 135°) and lowest temperature (-262. 8K @ 0 - 270°) are recorded near the 
heat exchanger entrance (x = 0 m.). 
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2-space model 


3-space model 



Figure 9(e-f). Temperature Profiles at Stationary Points. 

(Operating Conds.: 201.7 RPM, 1.008 MPa, T wall = 294 K, #tspc = 480, Grid size = 147 x 46 ). 


The 2-D codes peak temperatures are greater than Sage’s by ~ 40 K (2-space) and ~ 29 K (3 -space). The Sage code 
on the other hand is greater than the 2-D codes’ minimum temperature by ~ 25 K (2-space) and ~ 23 K (3 -space). 
The presence of the regenerator does appear to decrease the maximum temperature recorded in the heat exchanger of 
the 2-space model by ~ 11 K, and increase the minimum temperature by ~ 2 K; it also shifts the maximum and 
minimum temperature values from near the end of the heat exchanger to near the heat exchanger entrance during a 
cycle. 

The cylinder space and heat exchanger space pressure profiles are illustrated in Figs. 9(g-j). The pressure 
profiles in the cylinder and heat exchanger spaces are relatively insensitive to whether stationary or moving spatial 
coordinates are used; they are fairly symmetric over the crank angle range 0° < 0 < 360 0 with peak pressure values 
predictably occurring very close to the end of the compression phase of the cycle. In both cylinder and heat 
exchanger domains, Fluent and Sage pressure profiles are in excellent agreement over all crank angles for the 3- 
space model with maximum pressure 1.55E+06 N/m 2 occurring at 0 « 175°. For the 2-space model, CFD-ACE+ and 
Sage pressure profiles are not in agreement over all crank angles in the cylinder and heat exchanger domains. In the 
cylinder domain, the lack of correspondence is over the crank angle range 90° <<9<210° and in the heat 
exchanger domain, there is a significant difference between CFD-ACE+ and Sage pressure predictions. 
Correspondence is observed only at two crank angles (0^81° and 0 «261°). At maximum compression (0 « 175°), 
CFD-ACE+’s peak pressure prediction exceeds Sage’s by about 60.0 kPa with about 1.8° crank angle lead in the 
cylinder and by about 420.0 kPa with about 3.3° crank angle lead in the heat exchanger. It is also observed that 
whereas peak pressure values predicted by Fluent in the 3-space domains (cylinder and heat exchanger) are less than 
those predicted by CFD-ACE+ in the corresponding 2-space domains, peak pressure values predicted by Sage in the 
cylinder and heat exchanger domains are higher for the 3 -space model. Note also that in going from the cylinder 
space to the heat exchanger space, Sage predicts a large pressure drop (-0.36E+06 N/m 2 ). CFD-ACE+/Fluent 
predicts no pressure drop. The argument can be made that the addition of the regenerator helps to improve the 
pressure profile correspondence with the Sage result in both the heat exchanger and cylinder spaces and to reduce 
the 2-D peak pressure values in both the cylinder and heat exchanger domains. 
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Figure 9(g-j). Pressure Profiles at Stationary and/or Moving Points. 

(Operating Conds.: 201.7 RPM, 1.008 MPa, T wa n = 294 K, #tspc = 480, Grid size = 147 x 46 ). 


The pressure profiles in the regenerator are illustrated in Figure 10. Pressure profiles at two stationary points are 
recorded by the Fluent code. Fluent pressure profiles are insensitive to point location and are in excellent agreement 
with the pressure profile recorded by Sage over the crank angle range 0° < 6 < 360°. The maximum pressure of 
about 1.5E+06 N/m 2 occurs at 0 « 175°. 

Regenerator temperature profiles over the crank angle range 0° < 6 < 360 0 predicted using Fluent and Sage are 
illustrated in Figure 11. The fluid and solid temperature profiles predicted by Sage are in phase, each with very little 
temperature variation - only ~ 0.8 K around ~ 339.8 K. On the other hand, the fluid and solid temperature profiles 
predicted by Fluent are out of phase. The fluid temperature profile fluctuates about 40K around ~ 320 K and the 
solid temperature profile fluctuates about 37 K around 318K. There is clearly poor agreement between Sage and 
Fluent predictions. The arbitrariness of some of the prediction parameters (e.g., density and thermal conductivity of 
the solid matrix) introduced into the Fluent code in order to achieve the thermal non-equilibrium condition may be 
responsible for the disagreement between Fluent and Sage. As expected, the regenerator solid is seen to store heat 
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during the compression phase of the cycle (during which the fluid temperature is higher than the temperature of the 
solid) and to release heat to the fluid during the expansion phase of the cycle (during which the fluid temperature is 
lower than the temperature of the solid). 


Regenerator in 3-space domain 



Figure 10. Temperature and Pressure Profiles at Stationary Points. 

(Operating Conds.: 201.7 RPM, 1.008 MPa, T wa n = 294 K, #tspc = 480, Grid size = 147 x 46 ). 
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Figure 11. Temperature and Pressure Profiles at Stationary Points. 

(Operating Conds.: 201.7 RPM, 1.008 MPa, T wa n = 294 K, #tspc = 480, Grid size = 147 x 46 ). 
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1.3 Wall Heat Transfer Rate 

From the energy conservation principle expressed in differential form per unit time for a closed system 
undergoing a cyclic process: 

f5Q = c[<SW (47) 

The expansion and compression work of the piston 5W is a moving boundary work or pdV work where p is 
pressure and dV is differential volume. For non-quasi-equilibrium processes such as we have, the pressure at the 
inner face of the piston p f is used for p. Thus the equation relevant to our study is 

= Yfp- dv (48) 

The parameter T is the period of the cycle. In CFD-ACE+/Fluent, ^5 Q/t is calculated by first summing the heat 

transfer (Q) input data over the total number of time steps in a cycle and then dividing the sum by the period T and 
Jp.dV /t is computed using Simpson’s integration rule over the total number of time steps in a cycle viz: 



T 





(49) 


The superscripts n and n-1 imply values at the current and previous time steps respectively. The volume change 
(V n - V n l ) is calculated by multiplying the piston displacement during the given time step by the piston area. Sage 
output results provide mean values of the heat transfer rate and pressure which are equivalent to the integral 
parameters <j*6Q/T and (jp^v/T respectively. 

Recktenwald 27 modified Eq. (49) by introducing the work term weighting factor (3 to reduce the error introduced 
by the discretization of the pressure-volume work term viz.: 

j Vi dv= -V"-') (50) 

n=l 


The discretization is controlled by adjusting p to obtain cycle energy balance. As with the 2-space model, a p value 
of 1.49 was the optimum obtained for the 3 -space model. Table 6 shows work and heat transfer calculation results 
for a single cycle taken from mid 6 th cycle to mid 7 th cycle. 


Table 6. Sage and CFD-ACE+/Fluent Work and Heat Transfer Data 
(with the incorporation of the p weighting factor) 



Note that work and heat transfer calculations must be multiplied by two to get values for the entire domain since 
only one half of the domain is simulated. 
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The cylinder input heat rate corresponds to a rate of cooling of the cylinder walls. The heat exchanger output 
heat rate corresponds to the heating rate of the heat exchanger wall. As expected, the results indicate no energy is 
stored within the regenerator matrix during a cycle of operation. Generally, the heat stored in the regenerator matrix 
during one-half of the cycle is removed during the other half of the cycle. The net work input, w net in = Q out - Q in is 

confirmed by the results in Table 7. The piston action essentially pumps heat from the cylinder to the heat exchanger 
and regenerator. Thus with respect to the cylinder, the model acts like a “cooler” whereas with respect to the heat 
exchanger,, the model acts as a “heat pump” even though the wall temperature is kept constant at 294 K. 

With the addition of the regenerator, Sage reports ~ 23% increase in the cooling action in the cylinder, -11% 
reduction in the heat pumping action in the heat exchanger and - 25% reduction in the net cycle heat loss (or input 
work) is noted. On the other hand, Fluent, in relation to the CFD-ACE+ 2-space calculations, reports - 28% 
reduction in the cooling action in the cylinder (though the COP increases by - 12%), -33% reduction in the heat 
pumping action in the heat exchanger and - 36% reduction in the net cycle heat loss (or input work) is noted. Thus 
the addition of the regenerator leads to a reduction in the cylinder cooling action, heat pumping action and the net 
cycle heat loss (or work input) but an increase in COP. 

The heat addition-heat rejection process described above occurs in Stirling engines also. Heat is pumped from 
the expansion volume of a Stirling engine to the “appendix gap” (clearance volume between the displacer piston and 
the cylinder) and is lost to the walls of the clearance volume. This represents a net loss of heat to the work producing 
portion of the Stirling engine and is called an “appendix gap pumping loss”. 

Figures 12(a-d) show plots of wall heat transfer rate as a function of the crank angle. Cylinder wall, heat 
exchanger wall and total wall heat transfer rates are illustrated in Figs. 12(a,b) for the “2-space” model. The 
regenerator wall heat transfer rate is added in Figs. 12(c,d) for the “3-space” models. 



Figure 12(a-d). Wall Heat Transfer Vs. Crank Angle 
(Operating Conds.: 201.7 RPM, 1.008 MPa, T waU = 294 K, #tspc = 480, Grid size = 147 x 46) 
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Plots of the wall heat transfer rates exhibit oscillatory behavior over the range 0° < 0 < 360° in all the domains. 
It is observed that CFD-ACE+/Fluent and Sage signals are not in phase. In the 2-space model, the difference in the 
net cycle heat loss (or work input) between CFD-ACE+ (53.24 W) and Sage (29.95 W) predictions is ~ 43.7%. In 
the 3-space model, the difference in the net cycle heat loss (or work input) between Fluent (34.12 W) and Sage 
(22.37 W) predictions is ~ 34.4%. In each percent difference calculation, the 2-D code is used as the basis for 
accuracy. Note that these results are consistent with corresponding entries in Table 6. The incorporation of the 
regenerator tended to reduce the error between the 1-D and 2-D code predictions. 

2. CFD-ACE+/Fluent Results vs. Some Literature Results 
2.1 Surface Heat Flux and Temperature Difference 

Experimental and numerical heat exchanger surface heat flux and temperature difference between the gas 
temperature at the radial center of the heat exchanger and the heat exchanger wall at 1/16, 1/8, %, % of the heat 
exchanger length from the entrance and the “end” of the heat exchanger opposite the entrance are plotted against the 
crank angle starting at BDC. The surface heat flux data are illustrated in Figs. 13(a-d) and the corresponding 
temperature difference data are illustrated in Figs. 14(a-d). Komhauser’s experimental results 20 , Figs. 13(a) and 
14(a), Tew’s modified CAST code results 36 , Figs. 13(b) and 14(b), CFD-ACE+ results of the 2-space model 12 , Figs. 
13(c) and 14(c) and Fluent results of the present study of the 3-space model, Figs. 13(d) and 14(d) are presented and 
compared. Note that the sign of the heat transfer value generated by the modified CAST (Fig. 13(b)) is opposite to 
that of Komhauser’s (Fig. 13(a)). This is because in the numerical codes (CAST, CFD-ACE+ and Fluent), heat 
transfer is defined to be positive for heat flow from the wall to the gas (opposite to Komhauser’s definition). The 
signs of the heat transfer values generated by CFD-ACE+ (Fig. 13(c)) and Fluent results (Fig. 13(d)) were reversed 
in order to have the vertical coordinates appear like Kornhauser’s plot. The experimental (Kornhauser’s) and 
numerical (CAST, CFD-ACE+ and Fluent) results for the heat flux and temperature difference are generally in good 
qualitative agreement but with some notable differences. The variations in the numerical results appear smoother 
than Komhauser’s results. The magnitudes of the heat flux and temperature difference are also somewhat different. 
At maximum compression the highest values for the heat flux and temperature difference are observed at the 1/16 
and 1/8 location respectively for the experimental results and at the entrance of the heat exchanger for the numerical 
results. The heat flux and temperature difference values are smallest at the 1/8 position (Komhauser’s and Tew’s), at 
the end of the heat exchanger (2-space model) and at the % position (3 -space model). During the expansion phase 
(180° <6 < 360 0 ) when flow is from the heat exchanger to the cylinder, profiles of the heat flux and temperature 
difference generated by the CAST code are most insensitive to location along the heat exchanger surface. 
Komhauser’s results show significant profile differentiation for both the heat flux and temperature difference in this 
crank angle range. The minimum values are insensitive to position along the heat exchanger surface for the modified 
CAST code. There are also noticeable differences in the magnitudes of the heat flux and temperature difference 
reported by the codes at the heat exchanger entrance and end. 



Figure 13(a). Heat flux vs. Crank angle at various positions along 
heat exchanger surface relative to entrance to cylinder. 
(Komhauser’s Experimental Data 3 : Run #10271539, 201.7 RPM, 
1.008 MPa (arithmetic mean pressure), T wa u = 294 K). 



CRANK ANCLE FROM BC. DEC 

Figure 14(a). Temp, difference (T cent er- T wa ii) vs. Crank 
angle at various positions along heat exchanger surface relative to 
entrance to cylinder. (Komhauser’s Exptal Data 3 : Run #10271539, 
201.7 RPM, 1.008 MPa (arithmetic mean pressure), T wa n = 294 K). 
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Crank Angle from BDC, degree 

Figure 13(b). Heat flux vs. Crank angle at various positions 
along heat exchanger surface relative to entrance to 
cylinder. (Tew’s Modified Cast Code Calculations 70 *: 34 x 
20 srids. 120 tsnc.T 



Figure 13(c). Heat flux vs. Crank angle at various positions 
along heat exchanger surface relative to entrance to 
cylinder. (CFD-ACE+ Code Calculations 19 : 34 x 20 grids, 
120 tsnc.i 


Crank Angle from BDC, degree 

Figure 14(b). Temp, difference (T center - T wall ) vs. Crank 
angle at various positions along heat exchanger surface 
relative to entrance to cylinder. (Tew’s Modified Cast Code 
Calculations 70 *: 34 x 20 grids, 120 tsocT. 



angle at various positions along heat exchanger surface 
relative to entrance to cylinder. (CFD-ACE+ Code 
Calculations 19 : 34 x 20 grids, 120 tsnc.) 



Figure 13(d). Heat flux vs. Crank angle at various positions 
along heat exchanger surface relative to entrance to cylinder. 
(Fluent Code Calculations: 34 x 20 grids, 120 tspc). 


Figure 14(d). Temp, difference (T center - T wall ) vs. Crank angle 
at various positions along heat exchanger surface relative to 
entrance to cylinder. (Fluent Code Calculations: 34 x 20 grids, 
120 tspc). 
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The effect of the regenerator is clearly noticed at the entrance of the heat exchanger where the maximum values 
of the heat flux and temperature difference are elevated. 

2.2 Temperature Contours 

Temperature contour plots at 90° before TDC where maximum piston velocity is toward the heat exchanger are 
illustrated in Figs. 15(a-e). Similar qualitative results are obtained in the heat exchanger and cylinder spaces using 
the CFD-ACE+ code (Figs. 15(a-c)), the modified CAST code (Fig. 15(d)) and the Fluent code (Fig. 15(e)). It is 
observed that whereas in the heat exchanger space the fluid temperature increases from the wall region towards the 
interior (Fig. 15(b-e)), in the cylinder space the opposite trend is observed with fluid temperature decreasing as one 
moves into the cylinder interior from the wall region (Fig. 15(a,d,e)). Also, the cooler fluid is at the entrance of the 
heat exchanger and the fluid becomes warmer as one moves towards the end of the heat exchanger (Fig.l5(b,d,e)). 
Notice that unlike in the 2-space and 3-space models (Figs. 15(a,e)), the outer cylinder wall is flush with the outer 
wall of the heat exchanger in the CAST model (Fig. 15(d)). Due to this difference in geometries, it is observed that 
the higher temperature contours are situated close to the upper cylinder region in the 2-space and 3 -space models 
and not at the entrance to the heat exchanger as in CAST. 



Figure 15(a). Temperature contours for cylinder and 
entrance of HXer. (CFD-ACE+, 147x51 grids, At = 
6.19E-04, (480 tspc); 90° before TDC; Op. Conds.: 
201.7 RPM. 1.008 MPa) 


Figure 15(b). Temperature contours for middle part of 
heat exchanger. (CFD-ACE+, 147x51 grids, At = 
6.19E-04, (480 tspc); 90° before TDC; Op. Conds.: 
201.7 RPM. 1.008 MPa) 




Figure 15(c). Zoom in on temperature contours for end 
of HXer. (CFD-ACE+, 147x51grids, At = 6.19E-04, 
(480 tspc); 90° before TDC; Op. Conds.: 201.7 RPM, 
1.008 MPa) 


Figure 15(d). Temperature contours for entire domain 
(Modified CAST Run #13 with 82x20 grids, 960 tspc. 
Simulation of Komhauser’s [3] Experimental Run 
#10271539. Maximum Piston Velocity toward HXer 
90° before TDC.). 
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Figure : 15(e) Temperature contours of 3-space domain 
(Fluent, 147x46 grids A t = 6. 19E-04, (480 tspc); 90° befor TDC; Op. Conds.: 201.7 RPM, 1.008 MPa) 

As in the heat exchanger space, the fluid temperature in the regenerator increases from the wall region towards 
the interior. 


2.3 Velocity Vectors 

Figures 16(a-e) show laminar velocity vector field plots in the 2-space model (Figs. 16(a-d)) and 3-space model 
(Fig. 16(e)). The CAST code (Fig. 16(b)), the CFD-ACE+ code (Fig. 16(c)) and the Fluent code (Fig. 16(e)) show 
the fluid accelerating around the comer of the inner wall of the annulus to enter the heat exchanger. The velocity is 
primarily radial along the cylinder head approaching the heat exchanger entrance. Across the heat exchanger 
entrance, as the fluid prepares to turn the corner, the velocity changes from mostly radial at the lower “comer” to 
mostly axial at the upper wall for the CAST model. For the CFD-ACE+ and Fluent models, because of the slight 
difference in geometries, the transition from radial to axial flow takes place at the lower and upper comers of the 
heat exchanger entrance. Mass conservation requires a substantial increase in axial velocity near the outer wall for 
the CAST model (Fig. 16(b)) or in the middle, away from the walls for the CFD-ACE+ (Fig. 16(c)) and Fluent (Fig. 
16(e)) models, as the fluid enters the exchanger. Then the fluid begins to redistribute across the annulus and the axial 
velocity near the outer wall (CAST model) or near the inner and outer walls (CFD-ACE+ and Fluent models) 
decreases. 

As with the CAST results (Fig. 16(a)), the CFD-ACE+ results (Fig. 16(d)) show that the velocity tends to zero 
at the dead end of the heat exchanger in the 2-space model and Fluent results (Fig. 16(f)) show a reduction in the 
velocity field as the fluid moves toward the regenerator entrance. 
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Figure 16(a). Velocity field plot for entire domain 
(Modified CAST Run #13 with 82x20 grids, 960 tspc. 
Simulation of Komhauser’s [3] Experimental Run 
#10271539. Maximum piston velocity toward HXer.). 



Figure 16(c). Zoom in on velocity field plot at HXer 
entrance. (CFD-ACE+, Grid = 147x51, At = 6.19E-04, 
( 480 tspc); Maximum piston velocity toward the 90° 
before TDC; Operating Conditions: 201.7 RPM, 1.008 
MPa) 


Figure 16(b). Zoom in on velocity field plot at HXer 
entrance (Modified CAST Run #13 with 82x20 grids, 
960 tspc. Simulation of Komhauser’s [3] Experimental 
Run #10271539. Maximum piston velocity toward 
HXer.). 



Figure 16(d). Zoom in on velocity field plot at HXer 
end. (CFD-ACE+, Grid = 147x51, At = 6.19E-04, ( 480 
tspc); Maximum piston velocity toward the 90° before 
TDC; Operating Conditions: 201.7 RPM, 1.008 MPa.) 
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Figure 16 (e) Zoom in on velocity field plot at HXer entrance. Fluent Grid = 147x56, At = 6.19E-04, 
(480 tspc); Maximum piston velocity toward the 90° before TDC; Operating Maximum piston velocity 
toward the 90° before TDC; Operating Conditions: 201.7 RPM, 1.008 MPa). 
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Figure 16 (f) Zoom in on velocity field plot at HXer end. Fluent Grid = 147x56, 

At = 6.19E-04, (480 tspc); Maximum piston velocity toward the 90° before TDC; 

Operating Maximum piston velocity toward the 90° before TDC; Operating Conditions: 201.7 RPM, 1.008 MPa). 
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VI. Thermodynamic Loss Post Processing 

Equations (38-43) clearly show that local entropy production depends functionally on the local values of heat 
transfer rate, temperature, pressure, density, mass-specific entropy, velocity and viscous dissipation. Thus entropy 
generation can be considered a derived quantity that can be computed by post-processing experimental or numerical 
flow fields. Equation (44) shows that local entropy production is also a function of availability energy loss and 
temperature. 

1. Using Sage 

Sage entropy generation results (external and internal) for the 3 -space domain are calculated from Eq. (44) 
using availability loss result components obtained from Sage output file viz.: 


S_ = 


Availability Energy Loss AEfric + AEQw + AEQx ± | AEDiscr| 


(51) 


The parameters AEfric, AEQw, and AEQx are available energy losses due to flow friction, surface heat flow and 
axial heat flow respectively. |AEDiscr| is the absolute value of the discrepancy between the total available energy 
loss due to internal entropy generation (AEintemal = AEfric + AEQw + AEQx) and that due to external entropy 
generation, AEextemal. The “+” is used when it is assumed that AEextemal is greater than AEintemal and the “- 
“sign is used when the reverse is the case. These assumptions are arbitrary since the relative magnitudes of 
AEextemal and AEintemal cannot be determined apriori. A discrepancy of zero implies the total available energy 
loss due to internal entropy generation is equal to that due to external entropy generation. 


2. Using Fluent 

2.1 External Entropy Generation , S gen(ext) 


Closed System Analysis: 

Cylinder 

Piston 

Regenerator Heat Exchanger 


Axis of Symmetry 

Figure 17. 3-space Model Set-Up for External Entropy Generation Analysis (Closed System) 


The 3 -space model shown in Fig. 17 above taken as a whole constitutes a closed reciprocating system. When 
viewed as a closed system for external entropy generation analysis, Eq. (38) for the non-porous domain and Eq.(39) 
for the porous domain are thus reduced to: 


Non-porous domain: 


Porous domain: 


'gen,sys,cyl(ext) 

(non-porous) 


1 gen,sys, cyl(ext) 
porous 


/ 



dAdt 


dAdt 


> 0 


> 0 


(52) 


(53) 
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in which case only the results of surface heat transfer and temperature for both the fluid and solid phases are post- 
processed for external entropy generation. Eq. (52) is used to integrate the heat transfer along the isothermal 
bounding surfaces of the cylinder and heat exchanger and Eq. (53) is used to integrate the heat transfer along the 
adiabatic bounding surfaces of the regenerator. 

Open System Analysis: 

The 3 -space model can also be analyzed, by considering its separate components - regenerator, heat exchanger 
space and cylinder space - as open systems for which both entropy transfer components (heat and mass flow) of 
Eqs. (38) and (39) are now applicable. In order to simplify the analysis, the 3-space model was partitioned into five 
sub-domains - outer, mid and inner sub-domains in the cylinder space plus the heat exchanger and regenerator sub- 
domains as shown in Fig. 18 below. Fluent generated results of surface heat transfer, temperature, mass flow rate 
and mass specific entropy are post-processed for external entropy generation. 


Cylinder 



Outer 

S 

Piston 


Regenerator 

Heat Exchanger 

Mid 


Inner 


Axis of Symmetry 


Figure 18. 3 -space Model Set-Up for External Entropy Generation Analysis (Open System) 

The external entropy generation is calculated for each sub-domain with the sum of the results for the separate 
sub-domains yielding the total external entropy generation for the 3 -space model. The result for this open system 
analysis should equal the result for the closed system analysis. 

User Defined Functions 

User defined functions (UDFs) can be used to post-process any Fluent generated closed or open system data - 
surface heat transfer, temperature, velocity gradients, mass flow and mass specific entropy. UDFs are written in 
FORTRAN 90 and interfaced with the Fluent solver. The summation of surface, volume and time values, implied by 
the integrals in Eqs. (38-40 and 43) is facilitated with the use of do loops in the UDFs. In order to minimize program 
complexity, the UDFs are written to sum up the surface and volume integral values at each time step only and 
thereafter data is exported to an Excel spreadsheet for calculation of the cyclic time integral values. Carefully 
constructed algorithms are needed for implementation of Eqs. (38-40 and 43). These algorithms constitute the 
UDFs. 

For the closed system analysis, integrated surface heat transfer rates generated each time step are exported to an 
Excel spreadsheet for calculation of the cyclic time integral value. Division of the cyclic integral value by the 
constant surface temperature completes the calculation of the external entropy generation. 

2.2 Internal Entropy Generation, S gen ( int ) 

The sub-domains in the model may be further refined for calculation of internal entropy generation as shown in 
Figure 19, below, in order to obtain a better resolution of the losses. Internal entropy generation due to conductive 
heat flow and viscous dissipation is calculated for the non-porous region using Eq. (40) and internal entropy 
generation due to conductive heat flow, film heat transfer, viscous and inertial losses for the porous region using Eq. 
(43). Results of heat transfer, temperature and gradients of temperature, velocity and gradients of pressure are post- 
processed via UDFs written for each sub-domain in the 3-space model. 

With the entropy generation known, the availability energy loss can be calculated using Eq. (44). Focal and 
global distribution of entropy generation rates due to the various system irreversibilities can be evaluated via a 
histogram map. 
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Cylinder 



Figure 19. 3-space Model Set-Up for Internal Entropy Generation Analysis 


3. Entropy Generation and Availability Loss Results 

The results of the post-processing analysis are shown and discussed below. Figures 20 and 21 below illustrate 
the results of the external- and internal-entropy generations for each of the sub-domains in the 2-space and 3-space 
models. Figs. 20(a,c) show the functional dependence of entropy generation on the number of cycles of piston 
motion and Figs. 20(b,d) show histograms of the entropy generation values characterizing the contributions of each 
sub-domain to the external entropy generation. Figures 21(a-d) illustrate corresponding results for the internal 
entropy generation. 

Figs. 20(a,c) show the external entropy generation plots to be independent of the number of cycles of piston 
motion beyond the third cycle. The negative values for the external entropy generation in the cylinder sub-domains 
(Figs. 20(a-d)) should be interpreted as the cylinder acting as entropy sink. That is, these sub-domains appear to 
extract entropy from their surroundings (due to heat entering the cylinder). Because of the cylinder acting as entropy 
sink, the heat exchanger’s contribution exceeds the total external entropy generation by about 81% in the 2-space 
model (Fig. 20(b)) and by about 181% in the 3-space model (Fig. 20(d)). Figs. 21(a,c) show the internal-entropy 
generation plots to be independent of the number of cycles of piston motion one cycle earlier than does the external 
entropy generation plots. The heat exchanger’s contribution to the total internal- entropy generation falls short of the 
total by about 14% and 28% in the 2-space (Fig. 21(b)) and 3-space (Fig. 21(d)) models respectively. In this case, 
the cylinder sub domains contribute positively to the total entropy generation. 

Of the three cylinder’s sub-domains, the outer cylinder’s contribution to the entropy generation is most 
significant in the 2-space model (external or internal) and 3-space model (internal). The inner cylinder’s contribution 
to the external entropy generation is most significant in the 3-space model. The mid-cylinder’s and regenerator’s 
contributions to the external and internal entropy generations are minimal or negligible where relevant. When taken 
as a whole, the cylinder’s and regenerator’s contributions to both external and internal entropy generations are 
minimal compared to the heat exchanger’s contribution. 
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External Entropy Generation from 2-space Sub-domains (CFD-ACE+) 


Inner cyl 





Sub-Domain 


Figure 20(a). 2-space External Entropy generation vs. Cycle No. 
(201 .7 RPM, 1 .008 MPa., T wa u = 294 K ; Grid size = 147 x 46, #tspc - 480). 


Figure 20(b). 2-space External Entropy generation vs. Region 
(201.7 RPM, 1.008 MPa., T wa ii = 294 K ; Grid size - 147 x 46, #tspc - 480). 
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Figure 20(c). 3-space External Entropy generation vs. Cycle No. 
(201.7 RPM, 1.008 MPa., T wall = 294 K ; Grid size = 147 x 46, #tspc = 480). 


Figure 20(d). 3-space External Entropy generation vs. Region 
(201.7 RPM, 1.008 MPa., T waU = 294 K ; Grid size = 147 x 46, #tspc - 480). 
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Internal Entropy Generation in 2-space Sub-domains (CFD-ACE+) 
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Figure 21(a). 2-space Internal Entropy generation vs. Cycle No. 
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Figure 21(b). 2-space Internal Entropy generation vs. Region 


(201.7 RPM, 1.008 MPa., T wall = 294 K ; Grid size = 147 x 46, #tspc = 480). (201.7 RPM, 1.008 MPa., T wall = 294 K ; Grid size - 147 x 46, #tspc - 480). 
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Figure 21(c). 3-space Internal Entropy generation vs. Cycle No. 
(201.7 RPM, 1.008 MPa., T waU = 294 K ; Grid size - 147 x 46, #tspc - 480). 


Figure 21(d). 3-space Internal Entropy generation vs. Region 
(201.7 RPM, 1.008 MPa., T waU = 294 K ; Grid size = 147 x 46, #tspc = 480). 


A clear inference from these plots is that the major sub-domain contribution to entropy generation is from the 
heat exchanger. The heat exchanger is observed to contribute more to the total external entropy generation than it 
does to the total internal entropy generation. 

Figure 22 shows the distribution of entropy generation inside the 2-space and 3 -space models using Fluent. The 
impact of the entropy generation due to conductive heat transfer, fluid friction and mass transfer on the efficiency of 
the models is quantified. Under the specified condition (201.7 RPM, 1.008 MPa., T wa n = 294), relatively higher 
losses are experienced in the heat exchanger than in other sub-components of each model. This is followed by losses 
in the outer cylinder and inner cylinder. Contributions from the mid-cylinder are very negligible. The heat 
exchanger’s conductive heat transfer has the most impact on the model efficiency followed by outer and inner 
cylinder conductive heat transfer. Contributions to the entropy generation due to viscous dissipation, mass transfer 
and conductive heat transfer are very negligible. Please note that in Fig. 22(d), under the “Total Volume” column, 
“Viscous” has been mislabeled “Conductive” and vice versa. The relatively high losses in the heat exchanger 
support the extensive effort in the Stirling community to optimize the regenerator which is a form of heat exchanger. 
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In compact-heat-exchanger passages for example, improvements can be made in the constructive details of the 
channels (channel shape and aspect ratio, curvature of the return channels, etc.), the actual temperature differences 
between the cold and hot fluids, the surface finishing, and the type of materials used. 
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Figures 22. External and Internal Entropy Distribution in the 2-space and 3-space Models 
(201.7 RPM, 1.008 MPa., T wall = 294 K ; Grid size = 147 x 51, #tspc = 480, at Opt. Cycle = 6). 

Tables 7 and 8 illustrate entropy generation and availability loss results obtained for the 2-space and 3-space models. 
For the 2- space model, all CFD-ACE+ results, with the exception of the external entropy generation and external 
availability energy loss results for the cylinder, are greater than corresponding Sage results. For the 3 -space model, 
all Fluent results, with the exception of the external entropy generation and external availability energy loss results 
for the heat exchanger are less than corresponding Sage results. Calculation of the percentage differences between 
internal and external results of entropy and availability energy loss are based on the assumption that external results 
are more accurate since the calculation of the integral heat transfer rates are more dependable. With Sage, it does not 
seem to matter what the basis for calculation of the percentage differences is. Calculation of the percentage 
differences between Sage and the multi-dimensional (multi-D) codes (CFD-ACE+ and Fluent) are based on the 
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assumption that results generated by multi-D codes are more accurate since they are more suited to handle multi-D 
flow situations. Availability energy losses in the heat exchanger are greater than in the cylinder. The least losses are 
reported in the regenerator, as expected. The numerical errors reflected by the discrepancies between the internal and 
external entropy generation (or availability energy losses) indicate that Sage does a better job than the multi-D codes 
in satisfying the entropy generation accounting principle which require that the two methods of calculating for the 
entropy generation should give the same answer. 


Table 7: Entropy Generation and Availability Loss Results (Sage vs. CFD-ACE+) 


2-space (Cylinder + Heat Exchanger) 

Sub-Domain 

Sage 

(T wall = 294 K) 

CFD-ACE+: V2004 (Alpha version) 
(201.7RPM, 1.008 MPa., T wall =294K,) 
(Grid size = 147X51, #tspc = 480, Opt. Cycle = 6) 

Internal 

External 

Internal 

External 

Entropy Generation (KW/K) 

Heat 

Exchanger 

0.0000706900 

0.0000653500 

0.00013700 

0.000327000 

Cylinder 

0.00002245 

0.00002198 

0.000023000 

-0.000146000 

Total 

0.0000931400 

0.00008733 

0.00016000 

0.00018100 

Available Energy Loss (KW) | 

Heat 

Exchanger 

0.02120649 

0.01960549 

0.0402950 

0.096142 

Cylinder 

0.00673410 

0.00659400 

0.006793 

-0.042866 

Total 

0.02794059 

0.02619949 

0.047088 

0.053276 

AEDiscr 

0.00174110 

0.0061880 

Sub-Domain 

%Difference| between Internal and External Entropy Generation (or AE Loss) 

Heat 

Exchanger 

8.0 

58.0 

Cylinder 

2.0 

116.0 

Total 

7.0 

12.0 


|%Difference| between Sage and CFD-ACE+ Results 

Internal 

External 


Entropy 

AE Loss 

Entropy 

AE Loss 

Heat Exch. 

48.4 

47.4 

80.0 

79.6 

Cylinder 

2.4 

0.9 

115.1 

115.4 

Total 

41.8 

40.7 

51.8 

50.8 
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Table 8: Entropy Generation and Availability Loss Results (Sage vs. Fluent) 


3 -space (Cylinder + Heat Exch. + Regenerator) 

Sub-Domain 

Sage 

(T wa n = 294 K) 

Fluent 

(201.7RPM, 1.008 MPa., T wall =294K,) 
(Grid size = 147X51, #tspc = 480, Opt. Cycle = 6) 

Internal 

External 

Internal 

External 

Entropy Generation (KW/K) 

Regenerator 

0.0000067568 

0.0000064611 

0.0000000192 

0.0000012722 

Heat 

Exchanger 

0.0000352724 

0.0000313046 

0.00002645 

0.0001252 

Cylinder 

0.00002669 

0.00002270 

0.000010281 

-0.00008195 

Total 

0.000068718 

0.00006047 

0.00003675 

0.0000445222 

Available Energy Loss (KW) 

Regenerator 

0.002251 

0.00215253 

0.0000051835 

0.000343489 

Heat 

Exchanger 

0.01012557 

0.00898657 

0.0067493 

0.031940386 

Cylinder 

0.0076174 

0.0064784 

0.00262335 

-0.02020910121 

Total 

0.01999397 

0.0176175 

0.00937783 

0.01207477 

AEDiscr 

0.00237647 

- 0.00269694 

Sub-Domain 

%Difference| between Internal and External Entropy Generation (or AE Loss) 

Regenerator 

4.6 

98.5 

Heat 

Exchanger 

12.7 

78.9 

Cylinder 

17.6 

113 

Total 

13.5 

22.3 


|%Difference| between Sage and Fluent Results | 

Internal 

External 


Entropy 

AE Loss 

Entropy 

AE Loss 

Regenerator 

99.7 

99.7 

80.3 

84.0 

Heat Exch. 

25.0 

33.3 

75.0 

78.8 

Cylinder 

61.5 

65.5 

27.6 

67.1 

Total 

46.5 

89.4 

24.7 

35.4 


The errors for internal entropy and AE loss calculations between Sage and the multi-D codes are noted in 
general to be smaller than the errors for the corresponding external calculations. The exceptions are the errors 
reported in the regenerator. The great disparities in the percentage differences between Sage and the multi-D codes 
results of entropy and AE loss may possibly be due to the inability of the Sage 1-D code to accurately account for 
temperature and velocity gradients caused by flow separations from walls where there are changes in flow area as 
for example, between the cylinder and the heat exchanger in the 2-space and 3-space models. Sage 1-D flow is 
assumed to immediately adjust to the wall boundaries through all geometrical changes in area and can only 
approximately account for the effect of flow separations, vortex formation, etc. where there are changes in flow area. 
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VII. Conclusions and Suggestions for Future Work 

The following have been provided: an overview of the Stirling engine, a discussion of the computer models of 
the modified MIT 2-space test rig leading to the 3 -space model which includes a regenerator, the effect of the 
regenerator via comparison of the Sage 2-space and 3 -space modeling, theoretical development of the 
thermodynamic loss models, discussion of the numerical simulation results and post-processing of the numerical 
simulation results to obtain entropy generation/availability energy loss results for the 2-space and 3 -space models. 

The following conclusions can be drawn from the results of this study which sought to evaluate the effect of 
adding a regenerator to the MIT 2-space test rig and to characterize the irreversibilities related to heat transfer, mass 
flow and viscous friction occurring in the modified MIT test rig via entropy generation: 

(1) The inclusion of the regenerator: 

• reduces the minimum and maximum pressure and work-input values in the cylinder with the exception of 
the increase in Sage maximum pressure; 

• decreases the maximum temperature recorded in the heat exchanger of the 2-space model by - 11 K, 
increases the minimum temperature by - 2 K and shifts the maximum and minimum temperature values 
from near the end of the heat exchanger to near the heat exchanger entrance during a cycle; 

• leads, on the average, to - 22% reduction in the heat pumping action in the heat exchanger and -31% 
reduction in the net cycle heat loss (or input work) and whereas Fluent reports -28% reduction in the 
cooling action in the cylinder, Sage reports - 23% increase in same (both codes report increase in COP); 

• reduces the error between the 1-D (Sage) and 2-D (Fluent) code predictions; 

• elevates the cylinder space temperatures at stationary points especially at points close to the midpoint of the 
cylinder clearance volume (x = 0.001374 m) and piston top center position (x = 0.002747 m); 

• elevates the maximum values of the heat flux and temperature difference at the entrance of the heat 
exchanger; and 

• improves the pressure profile correspondence with Sage result in both the heat exchanger and cylinder 
spaces and reduces the 2-D peak pressure values in both the cylinder and heat exchanger domains. 

(2) The heat exchanger provides the major sub-domain contribution to entropy generation. Thus the heat exchanger 
heat transfer has the most impact on the 3 -space model efficiency followed by outer and inner cylinder heat 
transfer. Viscous dissipation throughout the entire 3 -space domain and mid-cylinder mass transfer and 
conductive heat transfer contribute minimally to the model’s efficiency. 

(3) In both the 2-space and 3-space models the two methods of accounting for the entropy generation (external and 
internal) appear insensitive to the number of cycles of piston motion beyond the third cycle. 

(4) The multi-D codes do not satisfy the 2 nd . Law accounting principle as well as Sage. These may be due to 
numerical losses. 

Clearly, there is room for further studies to extend the boundaries of this study and explore and improve upon 
the difficulties encountered in this study. For example, other regenerator parameter optimization studies such as wire 
diameter, annulus diameter (inner and outer) and porosity should be explored with the objective being the 
maximization of the cylinder-cooling effect. The difficulties encountered incorporating the loss models derived in 
Section IV in the Fluent code should also be revisited for possible resolution in order to provide a thermal non- 
equilibrium regenerator modeling capability in addition to the dual zone model provided by the Fluent code. 

Typically neglected and often viewed as superfluous, the second law of thermodynamics remains an esoteric 
and mysterious subject 15 particularly in computational analysis of thermo-fluid systems. When properly applied, the 
second law of thermodynamics has proven to be a very powerful tool in the optimization of complex thermodynamic 
systems. Loss analysis using entropy-generation rates due to heat and fluid flow is a relatively new technique for 
assessing component performance. It offers a deep insight into the flow phenomena, allows a more exact calculation 
of losses than is possible with traditional means involving the application of loss correlations and provides an 
effective tool for improving performance. Designers will know the cumulative amount of all losses computed locally 
in the flow domain. Entropy generation maps can be produced, and designers can use them by scanning them to 
detect critical areas (locations in which entropy generation is higher than its integral average value over the entire 
flow field). By considering the local values of entropy generation rates due to thermal and viscous dissipation, 
designers can generate a thermodynamically better design by simply trying to avoid these critical areas or re- 
computing them after a design modification has been introduced to assess local and global effects of the design 
change. 


NAS A/CR— 2008-2 1 5480 


46 



Our understanding of loss mechanisms is far from complete. Although numerical predictions are valuable in 
predicting the heat transfer and flow structure, there are difficulties in predicting the loss accurately. This is due to 
errors in predicting the boundary layers, transition as well as due to false entropy generation due to numerical 
dissipation. This work provides a point of reference for incorporation of loss post-processors into Stirling engine 
numerical codes. The incorporation of a loss post-processor in Stirling engine numerical codes, it is believed, will 
facilitate the optimization of Stirling engine performance. 
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